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Abstract 

The Schrodinger equation dictates that the propagation of nearly free electrons through a weak 
periodic potential results in the opening of band gaps near points of the reciprocal lattice known as 
Brillouin zone boundaries However, in the case of massless Dirac fermions, it has been predicted 
that the chirality of the charge carriers prevents the opening of a band gap and instead new Dirac 
points appear in the electronic structure of the material 2,l3|- Graphene on hexagonal boron nitride 
(hBN) exhibits a rotation dependent Moire pattern J, [5J. In this letter, we show experimentally 
and theoretically that this Moire pattern acts as a weak periodic potential and thereby leads to 
the emergence of a new set of Dirac points at an energy determined by its wavelength. The 
new massless Dirac fermions generated at these superlattice Dirac points are characterized by a 
significantly reduced Fermi velocity. The local density of states near these Dirac cones exhibits 
hexagonal modulations indicating an anisotropic Fermi velocity. 
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Due to its hexagonal lattice structure with a diatomic unit cell, graphene has low-energy 



electronic properties that are governed by the mass 
of consequences, among them Klein tunneling 



ess Dirac equation [6] . This has a number 



lOj which prevents electrostatic confine- 



ment of charge carriers and inhibits the fabrication of standard semiconductor devices. This 



has motivatec 
tials 



say 



a number of recent theoretical investigations of graphene in periodic poten- 



15| . which explored ways of controlling the propagation of charge carriers by 



means of various superlattice potentials. On the analytical side, one-dimensional potentials 
render particle propagation anisotropic 0, [j, ll|, [l^] and generate new Dirac points, where 
the electron and hole bands meet, at energies ±7m;f|G|/2 given by the reciprocal superlattice 
vectors G 

aa 

, and at zero energy for potentials with strengths exceeding m hvp\G\ [12|. 



Numerical approaches have extended several of these results to the case of two-dimensional 
potentials y, [l4, [l^. Unlike for Schrodinger fermions, the periodic potentials generally 
induce new Dirac points but do not open gaps in graphene, owing to the chiral nature of 
the Dirac fermions. 

Recent STM topography experiments have reported well developed Moire patterns in 
graphene on crystalline substrates, which suggests that the latter generate effective periodic 



potentials 



Il7j . Of particular interest is hexagonal boron nitride (hBN), because 



it is an insulator which only couples weakly to graphene. Furthermore, grap 



exhibits the highest mobility ever reported 
strongly suppressed charge inhomogeneities 



:or graphene on any substrate [18 
4 



rene on hBN 
and has 



. Hexagonal boron nitride is a layered 
material whose planes have the same atomic structure as graphene, with a 1.8% longer 
lattice constant. The influence of the weak graphene-substrate interlayer coupling on the 
electronic transport and spectroscopic properties of graphene is not well understood. In 
particular, there is to date no theory for local electronic properties such as those probed in 
STM experiments. Below we show that periodic interlayer couplings generate a new Dirac 
point at an energy determined by the wavevector of the periodic potential. The presence of 
this new Dirac point is reflected in two dips in the density of states, symmetrically placed 
at E = ±hvp\G\/2 around the E = graphene Dirac point but generally of asymmetric 
strength. There is also a periodic modulation of the local density of states, indicating 
an anisotropic Fermi velocity, with the same period as the superlattice topographic Moire 
pattern. 

The fabrication procedure used for creating the graphene on hBN devices results in a ran- 



dom rotational orientation between the graphene and hBN lattices. This rotation between 
the lattices and the longer lattice constant for hBN leads to topographic Moire patterns. 
Given the lattice mismatch 5 between hBN and graphene, the relative rotation angle (f) 
between the two lattices uniquely determines the Moire wavelength A as 

v/2(l + 5)(l-cos(0)) + (5 2 ' 

where a is the graphene lattice constant. The relative rotation angle 9 of the Moire pattern 
with respect to the graphene lattice is given by 

/„\ sin(0) . . 

tan {&) = - rN Vy; t~tt . 2 

w (1 + 5) - cosO) v ; 

Figure d^b) plots the wavelength of the Moire pattern (black) and rotation angle (red) as a 
function of <ft. Due to the lattice mismatch, there is a Moire pattern for all orientations of 
graphene on hBN with a maximum possible length of about 14 nm. Figs. [1] (c)-(e) show STM 
topography images of Moire patterns for three different rotations of the graphene lattice. 

We explore how this Moire structure influences the local density of states (LDOS) 
p(r,E) = J2 n IV ; n( r )| 2 ^(-^ — E n ) in the graphene layer. This quantity dominates dl/dV 
STM measurements as long as the density of states of the STM tip and the tunneling rate 
from the tip to the sample are constant. We take a lattice Hamiltonian for a graphene 
monolayer deposited on top of a hBN monolayer, where the two layers are rotated with 
respect to one another, and the spatially dependent interlayer hopping t± is calculated from 
nearest neighbor and next nearest neighbor interlayer coupling {4 ] . We numerically calculate 



p(r,E) for this model using the Lanczos method 19J (see Supplementary Information). 

Dips in the calculated p(r,E) are clearly seen in Fig. EJ^a). The energy of these dips 
chane; function of the rotation angle <fi and hence the Moire wavelength. We have 

also observed the dips in the experimental dl/dV curves as shown in Fig. [2](b). The black 
curve is for a 9.0 nm Moire pattern and the energy of the dip is 0.28 eV from the Dirac 
point. The red curve is for a 13.4 nm Moire pattern and the energy of the dips decreases 
to 0.22 eV from the Dirac point. Both experimentally and theoretically, we found that 
the relative strength of the dips in the conduction and valence band are different, with the 
dip in the valence band being much deeper than the dip in the conduction band. In our 
numerical calculations, we identified that most of this asymmetry arises due to next nearest 
neighbor interlayer coupling, which effectively induces modulated hopping between different 




Figure 1. Graphene device schematic and STM Moire images, (a) Schematic of the measurement 
setup showing the STM tip and an optical microscope image of one of the measured samples, 
(b) Superlattice wavelength A (black) and rotation (red) as a function of the angle between the 
graphene and hBN lattices, (c)-(e) STM topography images showing (c) 2.4 nm, (d) 6.0 nm and 
(e) 11.5 nm Moire patterns. Typical imaging parameters were sample voltages between 0.3 V and 
0.5 V and tunnel currents between 100 pA and 150 pA. The scale bars in all images are 5 nm. 

graphene sublattices and breaks electron- hole symmetry (see Supplementary Information). 
Figure [2(c) plots |d 2 I/dV 2 | for the 9.0 nm Moire pattern as a function of gate voltage and 
sample voltage. We clearly see the Dirac point in this measurement crossing the Fermi 
energy near zero gate voltage. There is a second dip which moves parallel to it that is offset 
by -0.28 V in sample voltage. This dip is due to the superlattice periodic potential induced 
by the hBN and indicates the emergence of a new superlattice Dirac point. 

We have observed these dips in the LDOS for seven different Moire wavelengths. The 
energy of the dips from the Dirac point is plotted (red points) as a function of wavelength 
in Fig. [2(d). The solid black line plots the expected energy dependence E = Hvf\G\/2 = 
2ixhv-p / \f3\ assuming the linear band structure of graphene and Vp — 1-1 x 10 6 m/s. For the 
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necessary high resolution spectroscopy, our STM is limited to observing dips in an energy 
range of ~ ±1 V which restricts the Moire wavelengths to longer than 2 nm. 




Energy (eV) ^ Sample Voltage (V) 




Gate Voltage (V) Wavelength (nm) 



Figure 2. Density of states of graphene on hBN showing new superlattice Dirac point, (a) Theo- 
retical LDOS curves for three different rotation angles between graphene and hBN, red is <p = 0.5° 
(12.5 nm), blue is 4> = 1° (10.0 nm) and green is <j) = 2° (6.3 nm). The curves have been ver- 
tically offset for clarity, (b) Experimental dl/dV curves for two different Moire wavelengths, 9.0 
nm (black) and 13.4 nm (red). The dips in the dl/dV curves are marked by arrows, (c) |d 2 I/dV 2 | 
as a function of gate and sample voltage for the 9.0 nm Moire pattern showing the shift of the 
Dirac point and one of the dips, (d) Energy of the dips away from the Dirac point as a function 
of Moire wavelength. The red points are the experimentally measured values and the black line is 
the expected theoretical dependence. 

To better understand these dips, we focus on the low-energy regime and neglect intervalley 
scattering in graphene. This is justified by the energy range in our STM experiments and the 
long wavelength of the Moire potential. The interlayer hopping term between the graphene 
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and hBN layers reflects the same periodic structure as the Moire pattern |4|]. Therefore we 
model the influence of the hBN by an effective periodic potential with the same symmetry 
as the observed Moire pattern. We accordingly consider the single-valley Hamiltonian 

if= hv F k-a + V^2cos(G a x)l, (3) 

a 

where k = (k x , k y ), a is a vector of Pauli matrices and X is the identity matrix. The potential 
strength is estimated as V = 0.06 eV from numerical second-order perturbation theory, and 
the G Q are the reciprocal superlattice vectors corresponding to the periodic potential gener- 
ated by the hBN substrate. The reciprocal superlattice vector Gi = (47r/\/3A)(cos0, sin#) 
is determined by the relative rotation of the graphene and hBN lattices according to Eqs. (JT]) 
and (j2J). The other superlattice wavevectors are obtained by two rotations of 60°. Larger 
superlattice vectors are not included in our model, since the corresponding couplings are 



smaller by more than one order of magnitude {4]. 

The dips in p(r, E) are due to k — > — k processes induced by the periodic potential for 
values 2k = G a corresponding to one of the reciprocal superlattice vectors G a . Unlike 
for Schrodinger particles, the chirality of the Dirac fermions prevents such processes from 
opening a gap at the edges of the superlattice Brillouin zone, as long as the potential does 
not break sublattice symmetry (see Supplementary Information). 

The presence of the new superlattice Dirac points can be detected by examining the 
gate dependence of the LDOS. As seen in Fig. [2](c), the two Dirac points move in parallel 
with gate voltage. Figure [3]^a) plots dl/dV for a 13.4 nm Moire pattern over a larger range 
of gate voltage. The white lines show the energy of the two Dirac points as a function 
of gate voltage. When the superlattice Dirac point approaches the Fermi energy, the gate 
dependence of the Dirac points changes. This is plotted in Fig. EJ^b), which tracks the energy 
of the main Dirac point when the superlattice Dirac point crosses the Fermi energy at -53 
V. Both Dirac points move more quickly with gate voltage, indicating a reduced density of 
states at this energy. The gate dependence for the main Dirac point can be fit with the 



equation E D = hv^^aix(y g — V a )/3 where V g is the gate voltage, V is the offset voltage 
and a is the coupling to the gate. The factor of three comes from the fact that the periodic 
potential creates three superlattice Dirac points corresponding to the three reciprocal lattice 
vectors G a . From the fit, we extract the value of the Fermi velocity at the Dirac point to be 
0.94 ±0.02 x 10 6 m/s for both electrons and holes. At the new superlattice Dirac points, the 
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Fermi velocity is reduced to 0.64 ±0.03 x 10 6 m/s for the new electrons and 0.78 ±0.03 x 10 6 

1 

m/s for the new holes. The reduction agrees with recent theoretical predictions [15j and 
our numerical calculations. Evidence for the presence of the superlattice Dirac point can 
also be seen in the global conductivity as a function of gate voltage (see Supplementary 
Information). 
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Figure 3. Gate dependence of graphene density of states near superlattice Dirac point, (a) dl/dV 
as a function of sample and gate voltage showing both the Dirac point and the new superlattice 
Dirac point. The white lines mark the location of the Dirac point and superlattice Dirac point, 
(b) Shift of the Dirac point as a function of gate voltage when the superlattice Dirac point crosses 
the Fermi energy (black line). The solid red line is a theoretical fit for the shift of the Dirac point 
with the presence of the superlattice Dirac point. The dashed red line shows the expected shift 
without the superlattice Dirac point. The inset shows the shift of both Dirac points over a large 
gate voltage range as well as theoretical fits (red and blue solid lines). 



In addition to the formation of superlattice Dirac points, the presence of the periodic 
potential also leads to a spatial variation in the LDOS. Figure H] shows experimental and 
numerical images of the LDOS as a function of energy. Near the Dirac point, (Figs. H] (b) and 
(e)), we observe a nearly featureless density of states. This is in agreement with previous 
STM measurements which show a strong suppression of charge fluctuations in graphene on 
hBN (4,5). At higher energies, the presence of the Moire potential manifests itself as a local 
variation in the density of states. Figures H] (a) and (d) are taken at a lower energy than 
the superlattice Dirac point in the valence band. At this energy the hexagonal pattern of 
the potential is clearly visible. On the conduction band side, Figs. IU (c) and (f), the Moire 



pattern is once again visible but its contrast is inverted. The centers of the hexagons now 
correspond to points of increased density of states. These LDOS modulations indicate that 
the Fermi velocity is spatially dependent. Moreover, the inversion indicates that the relative 
strength of the potential locally determines whether electrons or holes have a higher Fermi 
velocity. We observe that the Moire pattern becomes much more visible in the local density 
of states at energies above the superlattice Dirac points. Previous measurements of the 
LDOS in graphene on hBN did not observe any variations due to the short Moire patterns 




Figure 4. Experimental and theoretical images of LDOS for long wavelength Moire pattern, (a)-(c) 
Experimental dl/dV maps for a 13.4 nm Moire pattern. The sample voltages are (a) -0.16 V, (b) 
0.17 V and (c) 0.44 V. The sample voltage in (b) is near the Dirac point since the gate voltage 
was 60 V while the other two maps are near the energy of the superlattice Dirac points, (d)-(f) 
Theoretical dl/dV maps for a 13.4 nm Moire pattern. The energies are (d) -0.3 eV, (e) 0.03 eV 
and (f) 0.3 eV. The energy in (e) corresponds to the Dirac point. The scale bars in all images are 
10 nm. 

Graphene on hBN devices are becoming widely used due to their improved mobility and 
reduced charged impurities. We have shown that lattice mismatch and relative rotation 
between the graphene and hBN leads to a periodic potential for graphene charge carriers. 
This potential creates a new Dirac point whose energy is determined by the wavelength 
of the potential. This superlattice Dirac point has the potential to control the transport 
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properties of electrons in graphene as it induces spatially modulated velocities for the charge 
carriers. Future work is necessary to exploit this periodic potential for the creation of novel 
graphene devices. 



METHODS 

Graphene on hBN devices were fabricated using two different methods. In the first 
method, mechanically exfoliated graphene was transferred to high quality single crystals of 
hBN which were mechanically exfoliated on a SiC"2 substrate (4, 18 1. In the second method, 
commercially available hBN (Momentive AC6004) was exfoliated on SiO"2 substrates and 
then CVD grown graphene was deposited over the hBN. Both types of devices gave similar 
results so we do not distinguish between the two types. After depositing the graphene on 
hBN, Cr/Au electrodes were written using electron beam lithography. The devices were 
annealed at 350°C for 2 hours in a mixture of Argon and Hydrogen and then at 300 °C for 1 
hour in air before being transferred to the UHV LT-STM for topographic and spectroscopic 
measurements. We have measured a total of 29 samples of which 7 had long enough Moire 
patterns for the observation of a super lattice Dirac point. 

Figure la shows a schematic diagram of the measurement setup used for imaging and 
spectroscopy of the graphene flakes. All the measurements were performed in UHV at 
a temperature of 4.5 K. dl/dV measurements were acquired by turning off the feedback 
circuit and adding a small (5-10 mV) ac voltage at 563 Hz to the sample voltage. The 
current was measured by lock-in detection. 
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SUPPLEMENTARY INFORMATION 



I. LOW ENERGY PERTURBATION THEORY 

Scattering at the Brillouin zone boundaries due to the periodic potential does not open 
a gap for chiral massless Dirac fermions as long as sublattice symmetry is preserved. This 
can be seen by performing a unitary transformation H — >■ H' = UlHUi on the Hamiltonian 
of Eq. (3), with U\ = exp[zA(r) ■ a]. When V/hvp\G\ < 1 the linear part of the potential 
can be gauged out of the Hamiltonian for hv F V ■ A = — V. For V/hvp\G\ > 1, the gauge 
transformation can no longer remove the linear part of the potential, and the potential 
creates new Dirac points but still does not open a gap. jl-^]. 

The magnitude of the dips in the local density of states (LDOS) is determined by the 
strength and the wavelength of the Moire periodic potential. We see this perturbatively by 
projecting the Hamiltonian of Eq. (3) onto pairs of eigenstates 

■0s,k = (1, sexp^k])^ exp[ikr]/v^O 

of Hq = hvpla ■ a, with the lattice area Q and the angle #k = aictan(k y /k x ). The resulting 
2x2 Hamiltonian reads 

/ shv F k/2 V Kkl \ 

Hred = \ , (4) 

^ (V k>kl )* shv F k'/2 J 

with V k ,k' = (V/4)(l + exp[z(6» k / - 6> k )]) ]C a (<5k,G a +k' + <^k,-G a +k')- We focus on momenta 
k = G a /2 + 5k, k' = — G a /2 + 5k. For this choice of pairs of eigenstates, it is straightforward 



to see that the eigenvalues of if re( j are 

0(5k 3 /G 3 ), with e± = hvp\G a /2 ± <5k|, and the component 5k± of 5k that is perpendicular 
to G a . In particular, there is no energy shift for 5k || G a and a maximal shift for 5k - G = 0. 
This perturbative result is consistent with the numerical and experimental facts reported 
here that (i) no gap is opened, neither in the numerically obtained p(E,r), nor in the 
experimentally obtained STM dl/dV, (ii) a reduction of (dip in) the density of states is 
observed at an energy corresponding to HvfG/2, and (iii) both the spectroscopic Moire 
pattern and the dips tend to disappear for periodic potentials with shorter wavelength 
(hence larger G), when the relative rotation between the hBN substrate and the graphene 
layer is larger. The reduction of the density of states occurs around three points determined 
by the three reciprocal lattice vectors. At this level, there is no s-dependence of the strength 
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of the dip, the latter is the same for s = 1 (E > 0) as for s — — 1 {E < 0). An s- dependence 
emerges once off-diagonal terms ~ V'a x,y are included in the Hamiltonian of Eq. (3). 



II. BILAYER LATTICE HAMILTONIAN FOR GRAPHENE ON HEXAGONAL 
BORON NITRIDE 

We consider a lattice Hamiltonian for a graphene monolayer on a single layer of hBN. 
For each layer, the Hamiltonian reads 

H a = ^2 \ e A a ) a i( a ) a i( a ) + ^B{a)b\(a)bi(a) - t a ^(aj(a)6j(a) + h.c.) , (5) 

where a\(a) [a (a)] and b\(a) [6(a)] are creation [destruction] operators on sublattice A and 
B, respectively, of the graphene (a = 1) or hBN (a = 2) honeycomb lattice, and 
indicates that the sum runs only over nearest neighbors. On the hBN lattice, we choose the 
Boron atoms to be on the A sublattice and the Nitrogen atoms to be on the B sublattice. 
The on-site energies are e^(l) = £b(1) = 0, e^(2) = 3.34 eV and es(2) = —1.4 eV, and the 
hopping integrals are t\ = 3.16 eV and £2 = 2.79 eV [l]] (the precise value of the latter is of 
little importance). 

We restrict the interlayer hopping potential to nearest-neighbor and next-nearest- 
neighbor hopping. The interlayer hopping is given by 

t'ij( m > n ) = 7±exp[-|r;(m) - r.,-(n)|/£] f i} , (6) 

with a characteristic function fy = 1 if site i of sublattice m on the graphene (hBN) sheet 
is nearest or next-nearest neighbor to site j of sublattice n on the hBN (graphene) sheet, 
and fij = otherwise. The parameters r j± = ,0.39 eV and £ = 0.032 nm are calibrated to 



fit the interlayer couplings in bilayer graphene 



6] . The two lattices are rotated with respect 



to one another by the angle 0, and for each site we determine its nearest and next nearest 
neighbor site on the other sheet, numerically evaluate the distance between the sites and 
finally the corresponding interlayer hopping. While 7j_ should in principle depend on the 
sublattice index in the hBN layer, we neglect this dependence here. The interlayer coupling 
Hamiltonian is then given by 

^ = -E [4-(A^)[4(l)%(2) +h.c] + t' ij (A,B)[al(l)b j (2) + h.c.} 

ij 

+tUB,A)[b\(l) aj (2) + h.c] +4(S,B)[6l(l)6 i (2) + h.c.}] . (7) 
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Second-order perturbation theory maps 7i± onto a periodic potential of hexagonal symmetry, 
similar to the one in Eq. (3), with a modulation amplitude V ~ 0.06 eV which we determined 
numerically via second order perturbation theory in the interlayer hopping. This energy is 
smaller than Hvp\G\/2 with the Moire superlattice vector G, regardless of the rotation angle 
between graphene and hBN sheets. 

The total Hamiltonian reads Hi + H2 + H±- We evaluate the LDOS p(r i5 E) on the 
graphene sheet using the Lanczos method, which allows to reach linear system sizes of 
L = 1000 or more sites Q. The obtained p(r, E) depends on a smearing parameter ( which, 
as long as the density of states of the STM tip and the tunneling rate from the tip to the 
sample do not depend on energy can be related to the strength of the tip-graphene coupling. 
Figure ISlT a) illustrates how reducing this coupling allows finer and finer structures in the 
LDOS to be explored. For this particular rotation of <fi = 0.3°, corresponding to a Moire 
pattern with A = 13.4 nm, and the set of parameters we just discussed, we see that the 
LDOS vanishes more or less linearly close to E — Ed = —hvp\G\/2. Our numerical data 
suggest a complete, linear vanishing of the LDOS there. The linear vanishing of the LDOS, 
gives an estimate for the Fermi velocity close to the three new superlattice Dirac points. We 
estimate that the new velocity for the holes is (0.73±0.08)t>F- The new velocity for electrons 
is smaller because of the increased slope and we estimate it to be (0.53 ± 0.05)vp. The new 
reduced Fermi velocities are in reasonable agreement with our experimentally determined 
values. 

A significant energy asymmetry emerges in that the expected dip in the density of states 
is stronger in the valence E < than in the conduction E > band. This asymmetry arises 
because of (i) the asymmetry in the on-site energies in the hBN layer, which effectively 
induces a second order potential that is stronger for negative than for positive energies 
and (ii) next-nearest neighbor interlayer hopping, which effectively induces a local periodic 
modulation of t\. This is illustrated in Fig. ISTT b). 

III. TRANSPORT MEASUREMENT 

We have also performed electrical transport measurements on the graphene device shown 
in Fig. 3. The conductivity as a function of gate voltage is shown in Fig. IS2l The gate 
voltage is plotted as the offset from the gate voltage at the Dirac point. We observe two 
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Figure SI. Calculated density of states as a function of coupling, (a) LDOS for a 1000x1000 
graphene lattice on a hBN sheet, with relative rotation <j) = 0.3°. The smearing parameter of the 
Lanczos algorithm is C = 0.09 (black curve), 0.04 (red), 0.018 (green), 0.008 (blue) and 0.0036 
(violet). The dashed lines give linear fits close to the Dirac and superlattice Dirac points, (b) 
LDOS for a 1000x1000 graphene lattice on a hBN sheet, with relative rotation <j) = 2.0°, for 
Q = 0.0036. The black curve has nearest-neighbor interlayer coupling, the red curve has nearest 
and next nearest neighbor interlayer hopping. 

locations of decreased conductivity which are located at approximately ±40 V from the 
Dirac point. This is the same separation in gate voltage as observed in the spectroscopy 
measurements, where the main Dirac point was separated from the superlattice Dirac point 
by 40 V. Therefore, we conclude that the dips in conductivity are due to the presence of the 
superlattice Dirac point. This gives further indirect evidence of the superlattice Dirac point. 
We also see a second dip in the conductivity near -50 V. Our STM topography measurements 
showed a region of the device with a second Moire pattern of about 10 nm. The superlattice 
Dirac point due to this Moire occurs at a higher energy and hence a larger gate voltage. 
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Figure S2. Conductivity as a function of gate voltage for the graphene device with a 13.4 nm Moire 
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